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Abstract 

We describe quantum-field-theoretical (QFT) techniques for mapping quan- 
tum problems onto c-number stochastic problems. This approach yields re- 
sults which are identical to phase-space techniques [C.W. Gardiner, Quantum 
Noise (1991)] when the latter result in a Fokker-Planck equation for a cor- 
responding pseudo-probability distribution. If phase-space techniques do not 
result in a Fokker-Planck equation and hence fail to produce a stochastic 
representation, the QFT techniques nevertheless yield stochastic difference 
equations in discretised time. 



1 



Typeset using REVT^ 



I. INTRODUCTION 



There is a well-known duality between Fokker-Planck equations (FPE) and stochastic 
differential equations (SDE), which goes back as far as Einstein's and Langevin's theories of 
Brownian motion (see Risken's book for a detailed discussion of FP equation and related 
issues). Langevin equations or, more generally, stochastic differential equations, have long 
been a successful computational tool in quantum stochastics [^J^ , allowing for the numerical 
stochastic integration of systems for which analytical solution would be, at best, extremely 
difficult. There are long-standing rules, especially in quantum optics, for beginning with a 
particular system Hamiltonian and mapping this onto stochastic equations of motion for the 
field variables. However, this process depends on there being an FP equation equivalent to 
the master equation for the Hamiltonian in question; that is the partial differential equation 
for the appropriate probability or pseudoprobability distribution must contain derivatives of 
no higher than second order. This is the content of Pawula's theorem It restricts the 

method to a certain class of problems which, although containing many interesting cases, 
is not exhaustive. A generalised FP equation with higher order derivatives has no mapping 
onto a stochastic differential equation. As descriptions in terms of FP equations proper are 
the exception rather than the rule, it is of clear interest to extend present methods to allow 
for at least the numerical modelling of processes involving higher order noises. 

For purposes of numerical simulations, time may always be regarded as discretised, so 
the development of a stochastic difference equation approximating a generalised FP will 
generally be sufficient. Since Pawula's theorem only applies in the continuous time limit, 
this creates the loophole we need. What we will show here is that discretised stochastic 
equations may be devised for generalised FPEs, which, in a manner reminiscent of the 
positive-P representation of quantum optics, require a doubled phase-space. In agreement 
with Pawula's theorem, the method we use to develop these noises has no natural extension 
to the continuous time limit. This notwithstanding, any average corresponding to a physical 
quantity is expected to be correct in this limit (whereas the rest of the averages diverge). 
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The question whether a stochastic process in a certain generahsed mathematical sense can 
be defined corresponding to our methods is a subject for futher investigation. This is not, 
however, a problem for practical computer simulations, as these are necessarily performed 
on a discrete time grid. 

As in the well-known positive-P representation [0,0, nonphysical averages in our method 
determine the sampling noise. The fact that they are expected to diverge in the continuous 
time limit means that sampling noise grows as the time-grid spacing descreases. This is in 
contrast to the properties of the Wiener process, where sampling noise is independent of 
the grid spacing for a given sample size. In practice, however, this distinction between the 
Gaussian and higher-order noises is quantitative rather than qualitative, because for a given 
computing time the sample size necessarily decreases with the grid spacing. In either case 
the growth of sampling noise prevents the grid being made too fine. 

As a specific example, in this paper we consider a quantum oscillator with Kerr nonlinear- 
ity. It is well known that the Wigner function of this system obeys a generalised FP equation 
with 3rd order derivatives. Our goal is to show that the latter is (approximately) equivalent 
to a system of stochastic difference equations. This equivalence is not straighforward but 
rather bears a lot of similarity to the positive-P representation known in quantum optics 
PJ^. We shall therefore call these equations a positive-W representation of the nonlinear 
quantum oscillator. 

The positive-P representation originates with the single-time P-distribution, yet it allows 
one to calculate a much wider class of quantum averages 0,0], namely, averages of multi-time, 
time-normally ordered ||^ operator products. The positive-P representation may thus be re- 
garded as a constructive mapping of a quantum problem to a classical stochastic problem. 
Its characterictic property is that time-normal averages are mapped directly on classical 
averages. In a much similar manner, the positive-W equations, which originate with the 
single-time W-distribution, allow one to calculate multi-time, time- Wigner ordered operator 
averages. This new type of operator ordering is introduced in this paper; to the best of 
our knowledge it has never been discussed before. It generalises the single-time symmetric 
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operator ordering, in the same way as the time-normal ordering generahses the single-time 
normal ordering. The positive-W representation is hence a constructive mapping of a quan- 
tum problem onto a classical stochastic problem, with the characterictic property of relating 
time-Wigner averages directly to classical averages. 

The question then is if, rather than taking the usual route from q-number to c-number 
equations based on phase-space techniques the positive-P and the positive-W repre- 

sentations may be derived by directly linking q-number Heisenberg equations to c-number 
Langevin equations. In this paper, we show that such a derivation is indeed possible, by 
employing methods P,^] based on the techniques of quantum field theory. Our considera- 
tions closely follow the way in which Feynman diagram techniques were derived in textbooks 
dating back to the fifties and sixties 0. We derive Keldysh diagram series for the Kerr 
oscillator, then recast them as a Wyld-type series |]10[, otherwise termed causal series [|7|,p 



[see also Ref. |T2[)- (More precisely speaking, we derive generating expressions [0 for these 



series, avoiding diagram notation as such.) Causal series emerge as solutions to c-number 
stochastic problems ||TT|. On the other hand, given a causal series, it is easy to write a 
stochastic differential (or difference) equation for which this series is a solution |]TT|. This 



yields strikingly simple and powerful techniques for obtaining stochastic representations of 
quantum problems. (Cf., e.g., section [VB 1| where positive-P equations are derived for an 
optical parametric oscillator, the essence of the derivation being a two-line calculation.) Not 
the least important property of these techniques is that they can be formulated using simple 
recipes and then used without any reference to the advanced methods employed in their 
derivation. 

The paper is structured as follows. In section we reiterate the standard techniques 
of quantum field theory as applied to a nonlinear oscillator. We discuss in detail causal 
regularisation |Tl[] of the propagator, which is necessary to order to make our relations 



unambiguous. A functional (Hori's) form of Wick's theorem [jl3[ is introduced for the normal 
ordering and then generalised to the case of symmetric ordering. The result of the section 
is a pair of closed perturbative relations. One of them applies when the inital state of the 



field is characterised via a P-distribution; the other apphes if this state is characterised via 
a W-distribution . In section ^TT|, we investigate the causal structure of the perturbative 
relations derived in section |I|. The concept of causality is introduced via the retarded Green's 
function of the free Schrodinger equation. Following it, we are able to define the input and 
output of a quantum system. We then show that, physically, the input and output thus 
introduced correspond to a generalisation of Kubo's linear reaction approach to a full 



nonlinear quantum-stochastic response problem. In section we develop techniques (based 



on the Hubbard- Stratonovich transformation [|T5l) which allow for a constructive mapping 
of the quantum nonlinear response problem onto a classical nonlinear stochastic response 
problem. For the oscillator, this results in the positive-P and positive-W representations 
which we are seeking. We complete the section with results of computer simulations of the 
Kerr oscillator using these representations. Finally, in section we reformulate our results 
recipe-style, the way they should be applied in calculations. Despite all the tediousness of 
their derivation, we end up with two simple relations, which allow one to derive positive-P 
and positive-W representations for an arbitrary nonlinear quantum system. We illustrate 
their use by applying them to a number of systems common in quantum optics. 



II. QUANTUM FIELD THEORY OF THE KERR OSCILLATOR 

A. The model 

The techniques we introduce in this paper are applicable to any Hamiltonians which 
have a polynomial form in the field operators. We also assume that the Hamiltonian can 
be divided into a quadratic part, called the free Hamiltonian, and the remainder, termed 
the interaction Hamiltonian. This allows us to introduce, in the usual manner, Schrodinger, 
Heisenberg, and interaction-picture field operators. However, this requirement may always 
be satisfied by subtracting a suitable quadratic term from the Hamiltonian, and declaring 
the remainder as being the "interaction" part (cf. the way a Higgs-type phase transition in 
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an anharmonic oscillator was treated in Ref. |^). 

This generality notwithstanding, our techniques may be effectively demonstrated for a ID 
oscillator (as is usual in quantum field theory). We therefore consider a nonlinear quantum 
oscillator with the Hamiltonian, 

n = no + Hint = t^fls^s + J al^al, (1) 

using units such that h = 1. In Og and as are the usual pair of creation and annihilation 
operators with commutator [as,as] ~ ^- They play the role of Schrodinger-picture field 
operators for this illustrative system. The field operators in the interaction picture are 
simply 

a{t) = e-'^^'as, a\t) = e*'^*a^, (2) 
while Heisenberg picture operators will be denoted in a slanted font as a'it) and a(t). 



B. Time orderings of operators 

Time ordering of operators puts operators from right to left in the order of increasing 
time arguments, e.g., 

I at(t)a(t'), t > f, 
T+a{t')a\t) ^ J ^ ^ ^ ^' - ' (3) 

I a{t')si\t), t < t'. 

For equal times, the time ordering is specified as normal ordering (which places all creation 
operators on the left of annihilation operators). Further, the reverse time ordering, T_, 
places operators in the order of decreasing time arguments. Formally, it may be defined as 
a conjugate of T+: 



where P is a product of field operators. Then, e.g.. 



tUp')]'^ (4) 
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, aUt)a(t'), t < t', 
T_k{t')k\t) _ ^ ^ ^ ^ ^' - ' (5) 

a{t')ei\t), t > t'. 

For equal times, T_ also becomes normal ordering. Double time ordering is the combination 
of the T+ and T_-orderings, 

T_P_ (6) 
where P_ and P+ are operator products (by definition, T_ acts on P_ and T+ acts on P+). 

C. Non-stationary perturbation approach 

1. Perturbation expressions for time- ordered operator products 

Heisenberg operators are related to those in the interaction picture via the evolution 
operator, 

a{t) = U\t, -oo)a{t)U{t, -oo), (7) 

a}{t) = U\t, -oo)a){t)U{t, -oo), (8) 

which is a solution to the Schrodinger equation. 



dt 

where 



= nint{t)u{t, to), u{to, to) = a, (9) 



nUt) = ^a^'{t)a'{t). (10) 
The evolution operator may be written as a time-ordered operator expression (T-exponent) : 



-i fdt'HUt') 

Jto 



(11) 



U{t, to) = T+ exp 
It is unitary, and has the group property, 

U{t, t')U{t', t") = U{t, t"), t>t' > t". (12) 



In particular it follows that 



—oo) 



(13) 



Then, assuming that times are ordered in a certain way, t > t', we have: 



T+a\t)a{t') = U\oo, -oo)W(oo, t)a\t)U{t, t')a{t')U{t' , -oo). 



(14) 



It should be noted here that the combination 



W(oo, t)a\t)U{t, t')a{t')U{t', -oo) 



(15) 



is time-ordered. Finally, introducing the S-matrix, 



S = W(oo, — oo). 



(16) 



we have. 



T+a^(t)a(t') = S^T+a{t')a\t) exp 



IK 



+ 00 



r a r 



(17) 



There is a certain subtlety in the way in which relations such as ( |Tl| ) and ( p!7D should 
be understood. Quantities under ordering should be regarded as functions of d(t) and a"'"(t), 
whereas (0) may be used only after the time ordering has been completed. So calculating 
the time integral in (|T^ should follow the ordering rather than precede it. Were one to 
take the integral before the time ordering, it would prevent T+ from acting on the operators 
comprising the interaction Hamiltonian. 

2. Functional perturbation techniques 

For our purposes we need relations applying to the whole assemblage of time-ordered 
operator products. An arbitrary time-ordered operator product may be found by differenti- 
ating the following time-ordered operator exponent. 



T+ exp 



dt 



C{t)a^{t) + CHt)a{t)] ]=T+ exp (Ca^ + C^a) , 



(18) 



where ((t),(K't) is a pair of c-number functions. Then, e.g., 

T^a{t')a\t) = exp ((at + C^a) |^=^t=o, (19) 

and so on. Equation {\L7\} is readily generalised, resulting in 



r+ exp (Ca^" + C^a) =5%exp|y dt C{t)d\t) + C\t)d{t) - — a^\t)d\t) 



S^T+ exp [Ca^ + C^^ " j ■ (20) 



Relations for the inverse and double time orderings then read 



T_ exp (Ca^ + C^a) = T_ exp (^(d^ + d^^d^j S, (21) 

T_ exp (^C_a^ + Cia) T+ exp (c+a^ + (laj 



= T_ exp (^C-a^ + + ^ d^^d^j T+ exp (^C+d^ + C+a - ^ j • (22) 

The relation for the double time ordering (naturally) involves four arbitrary c-number func- 
tions, C±(^)5C±(^)- Note that, in this relation, the factors S and 5^ outside the orderings 
have cancelled each other. This results in a genuine double-time-ordered structure on the 
RHS of (0). 

In what follows, we will make wide use of the brief notation as in Eqs. (|lB), pO|) and 



(p^), which will always be signalled by the absence of time arguments of the fields. Apart 
from these cases, we will not omit time arguments where they apply, as in Eq. (0), so as to 
remove any ambuguity. (Note that this is also the reason why the Scrodinger-picture field 
operators were denoted as,ag, not a,a^) 



D. Wick's theorems 

1. Hori's form of Wick's theorem proper 

A common way of dealing with time-ordered expressions is by converting them to nor- 
mally ordered form. This is done following Wick's theorem, which states that "a time ordered 
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product of interaction-picture field operators equals the sum of all possible normally ordered 
operator products, obtained by replacing pairs of operators in the initial product by cor- 
responding contractions (including the term without contractions)". For the oscillator, we 
have only one nonzero contraction, namely. 



T+a{t')a\t)- : a\t)d{t') : 







T+d\t)a{t') 







(23) 



In the above, G{t) is a retarded Green's function of the free Schrodinger equation: 

d 



It then follows that 



T+d\t)d{t') 
T+d^{t)d{t')d{t") 
T+d\t)d\t')a{t")d{t"') 



^--uJ\G{t) = 6{t). 



d\t)d{t') +iG{t' -t), 

d^{t)d{t')a{t") + iG{t' - t)a{t") + iG{t" - t)d{t'), 
a\t)a\t')a{t'')a{t''') 

+ iG{t" - t)d\t')a{t"') + %G{t"' - t)a\t')a{t") + 
+ iG{t" - t')d\t)a{t"') + iG{t"' - t')a\t)a{t") 
- G{t" - t)G{t"' - t') - G{t" - t')G{t"' - t), 



(24) 



(25) 
(26) 



(27) 



and so on. 

As was noticed by Hori [|1^, the pattern of products with contractions is exactly that 
produced by a quadratic form of functional derivatives. 



(28) 



acting repeatedly on products of c-number functions a(t), a^{t). Applying A once produces 
terms with a single contraction: 



Aat(t)a(t') = iG{t' -t), 
Aa^{t)a{t')a{t") = iG{t' - t)a{t") + iG{t" - t)a{t'), 



(29) 
(30) 
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Aa\t)a^{t')a{t")a{t"') = iG{t" - t)a\t')a{t"') + iG{t"' - t)a\t')a{t") + 

iG{t" - t')a\t)a{t"') + iG{t"' - t')a\t)a{t"), (31) 

etc. Terms with n contractions are produced by A"/n!. As an example, consider the case 
where n = 2: 



A2 



— a\t)a\t')a{t")a{t"') = - G{t" - t)G{t"' - t') - G{t" - t')G{t"' - t). 



(32) 



The whole assemblage of terms required by Wick's theorem (including those without con- 
tractions) is thus found by applying a differential operator, 



. A A' 

1 + A + - + .. 



(33) 



This allows one to write Wick's theorem in a compact form as (with P being an arbitrary 
operator product) 



-'^ {P\a{t)-*a(t),a 



t(t)-.at(i) 



la(t)^a{t),at(t)-^at(t) 



(34) 



2. Causal regularisation 

Wick's theorem requires that no contractions should occur between operators with equal 
time arguments. (This is simply because the time ordering was specified for equal times as 
normal ordering.) However, applying ( P^ results in ambiguities, e.g., 

Aa^{t)a{t) = iG{0) = ie{0) (undefined). (35) 

A convenient way around this problem is a causal regularisation of G{t). To this end, we 
shall always assume that G{t) is somehow smoothed while still preserving its causal nature: 
G{t) = for t < 0. For a continuous function this also means that G{0) = 0. With this 
specification (and implying a final limit of unsmoothed G), Eq. ( ^If ) becomes fully equivalent 
to Wick's theorem. 



11 



3. Wick's theorem for double time ordering 



It may be checked that proof of Wick's theorem ^ is based only on the hnear ordering of 
the time axis; consequently Wick's theorem may be generalised to operators defined formally 
on any linearly ordered set. This clearly applies to the double time ordering, Eq. (P), which 
may alternatively be introduced as an ordering on the so-called C-contour ^j. The C-contour 
(see Fig. |ID first travels from t = — oo to t = +00 (direct branch) and then back to t = —00 
(reverse branch). Then, 



I a — * a-f , a I — *a 



(36) 



with the subscripts '+' and '— ' specifying operators as being defined on the direct and re- 
verse branches of the C-contour. (In other words, operators acquire an additional argument; 
a C-contour index. This argument is useful purely for labelling purposes and should be dis- 
regarded after the ordering has been performed.) For the Tc ordering, operator contraction 
becomes a matrix in respect of the C-countour indices, (a,/5 = -|-, — ) 



iGap{t' - t) = (0 Tcao,{t')a\{t) 







(37) 



and the functional (Hori's) form of Wick's theorem generalised to the double time ordering 
is found to be: 



T_P_ T,P, 



Ac p_ 



I a{t)^a_ (t),at (t)^a'_ (t) 



X P_ 



+ I a(t)^a+ (t),at {t)-^a}, (t) 



(38) 

Here, P_ and P+ are arbitrary operator products, a±(t), aj_(t) are four independent c- number 
functions, and Ac is a quadratic form of functional derivatives, 

52 



dtdt' 



5a+{t')5a\{t) 



+ iG__(t'-t)- 



52 



+ iG_+(t'-t)- 



6a.{t')6al{t) 8a-{t')6a\{t) 
where the kernels are the three nonzero components of Gap (shown schematically in Fig. 
In turn, these are conveniently expressed in terms of G{t): 



(39) 
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G++{t) = G{t), G__(t) = - G*{-t), G_+(t) = G{t) - G*{-t). (40) 

Note that Eqs. (|38D and ( ^O]) are formulated so as to eliminate the concept of the C-contour 
from our considerations from here on. 

Relations (^OD require a word of caution. For equal times, both the and the T_ 
orderings are specified as normal ordering, hence the no-contractions-between-the-same- 
time-operators caveat of Wick's theorem applies equally to operators under the T_ ordering. 

In (0), it is enforced through the equation relating G (t) to G(t) and the (implied) causal 

regular isat ion of G(t). At the same time, through the equation relating G—^-{t) to G{t), the 
regular isat ion also modifies G_+(t), "burning a hole" in it in the vicinity of t = 0. Both 
the causal regularisation and relations (^OD play central roles in our considerations, so we 
have to make sure that the stated modification of (^-^(t) does not lead to incorrect results. 
The fact that weird results may indeed follow demonstrates, e.g., a "proof" that ds and Sg 
commute. With regularisation, G—^{0) = 0, and we "obtain", 

asal = Tca-(t)at W = a\t)a{t) + 2G_+(0) = a^ag. (41) 

The flaw in this "proof" is that all quantities we deal with should be regarded as generalised 
functions (distributions) and not pointwise functions. This is especially true under regu- 
larisation. "Holes" in continuous functions which emerge due to regularisation should be 
simply ignored (smeared out). We would only expect problems associated with the "holes" 
if they were overlapping with sufficiently strong singularities, 5-functions or worse, whereas 
the worst type of singularity that we may expect to occur is a step-function. The "hole" in 
G_+(t) should thus be of no consequence. 

4- Reordering time-ordered operators symmetrically 

Wick's theorem also holds if the desired result is symmetric, or Wigner, operator ordering. 
(For a definition of the symmetric ordering see, e.g. [0.) The easiest way to demonstrate 
this is to derive an analog of Eq. ([381) . A relation between a normal and a symmetric 
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representation of a particular Scrodinger operator, most suitable for our purposes, is given 
in Ref. |T^. It can be readily adapted to operators in the interaction picture resulting in 

S ^ 6 



: P : = W 



la{t)^a(t),at(t)-+at(t) 



(42) 



Here, P is an arbitrary operator product, W denotes Wigner ordering, and 

s s s'^ 



(43) 



where the kernel Gw is defined as 

iGw{t - t') = d\t)d{t') - Wa{t')a\t) 







Wd{t')a^{t) 







(44) 



Note that the type of brief notation introduced by (|43|) accounts for asymmetry of the kernel 
Gw- We now make use of the following fancy form of the rule of product differentiation. 



$ {£\ $1 (V^) (^) = $ + ^) *1 M $2 (^.) 1.,,..^.. 



(45) 



where (/^(t), V'i(t), V2if) are c-number functions and ^{^), ^i{ip), $2(v') are functional of 
such functions. Combining (133) and (|42D and employing (03), we arrive at a relation. 



T_P_ T, P, = W 



e^c P_ 



I a(t)^a_ (t),at (t)^al (i) 



X P 



+ 1 a{t)^a+ [t) ,at (t)^a;'|_ (t) 



la_(t),a+(t)^a(t),al(t),a^(t)->at(t) f ' 



(46) 



where 
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(47) 



For nonzero time arguments, kernels comprising ( |77| ) may be expressed by the retarded 
Green's function, G{t). We note that 



Gap{t) + Gw{t) = G^p{t) - -[G{t) - G*{-t)] ^ G^p{t). 



(48) 



Then, employing (140|), we find. 



G^_,(t) = -G^_(t) = \ [Git) + G*(-t)] , G^^t) = -G^_(t) = \ [Git) - G\-t)] . (49) 
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Extending these relations to all times, we define, 

A?'=» E i^G- * • (50) 

The question then is if in (^H]) may be replaced by . The same arguments as above 
show that the "holes" in G^_^_{t) and G^_{t) are of no concern. This is not the case for 
G^^{t) and G^_{t). Indeed, recall that the and T_-ordering are both specified for equal 
times as normal ordering. Aq takes care of such same-time operator groups reordering them 
symmetrically, whereas misses them. That is, enforces the no-contractions-between- 
the-same-time-operators caveat of Wick's theorem, whereas this caveat does not apply if a 
double-time-ordered operator product is reordered symmetrically. We should thus either 
amend Wick's theorem introducing same-time contractions, or, which is more convenient 
and suitable for our purposes, redefine the time orderings. We then obtain 



e^^ ( P_ 



T^p_ r_f P. = w 



^ 1 a(t)^a+(t), at (t)^at(t) 



(51) 



la_(t),a+(t)^a(t),aL(t),a^(t)^at(t) f ' 

were and differ from T_ and T+ in that same-time operators are ordered symmet- 
rically rather than normally. 



E. Closed perturbative relations for quantum field averages 

1. The method of normal ordering 

For practical purposes, the quantities of interest are operator averages rather than the 
operators themselves. We therefore consider a characteristic functional of averages of double- 
time-ordered operator products: 

H (C_, C+, Cl Cl) = exp (C-at + (la) T+ exp (c+at + C^a) ) . (52) 

Angle brackets here define an averaging over the Heisenberg p-matrix of the quantum field, 
(or over the field's initial state, which is the same thing): 
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•) = Trp(---). (53) 



Although defined in q- number terms, in itself functional ( [5^ ) is a c-number object. Ap- 
plying the closed perturbative relation, Eq. (P^D, and then Wick's theorem, Eq. (pSD, allows 
one to express it as an average of a normally ordered operator expression. To eliminate q- 
numbers completely, we shall characterise the initial state of the field by the corresponding 
P-distribution, 

P(a) = ^Jd^V (e<^''^'y^'^^"-)-\^^\'l^'^ , (54) 
p = cfaP{a) \a) {a\ , (55) 



where \a) is a coherent state with the amplitude a: 

a\a) = a \a) . (56) 
For any normally ordered operator expression, : X :, we then have 

[: X :) = / £aP{a) la:X:a)= d^aP{a)X\a(t)^a{t),dHt)^a*{t), (57) 



where a{t) = ae is the coherent amplitude of the interaction-picture field operator, 
d{t) \a) = a(t) I a). Combining Eqs. (p^), ( p8| ) and ( |57D then yields: 

\a-it)=a+{t)=a{t),al{t)=a''{t)=a*{t)- \'^°) 



2. The method of Wigner ordering 



In order to derive an analog of Eq. (|58|) for symmetric ordering, consider, to start with, 
relation ( pT]) for the evolution operator. By derivation 0, it is a sum of T-ordered terms, 
found by expanding the exponent in a power series: 



U{t, to) = Texp 



-t l' dt'Hintit') 
'to 



n-t dt'nut')+ o 



to 



dt'dt"Tnut')'>iUt") + 



(59) 
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The T-ordering here apphes, strictly speaking, to the Hamiltonian operators regarded as 
entities rather than to the field operators. Redefining the time ordering for the field operators 
via Eq. (Rf) changes the integrands, because, e.g.. 



(60) 



This change, however, (i) is finite and (ii) affects only the submanifold of measure zero of 
the integration manifold. Hence it has no effect on the result of integration (this and further 
arguments below may become flawed for continuous-space problems where expressions like 
['^int(^)]^ — : ['^int(^)]^ '■ often contain inflnities). 

The key property of the T+-ordering is that it does not affect the Hamiltonian operators 
themselves, 

r+{?^i,t(t)}=Hint(t), (61) 

which in turn is due to the normal form of 'H\at{t\ It is then immediately clear that (pUf ) 
may be rewritten in terms of by using a symmetric form of the interaction: 



\a^\tra\t) = ^W{at^(t)a^(t) - 2a\t)ait) + ^} ^ <,(t). 



(62) 



so that 



When deriving (^2p we have used the property that 

W |a^^a^| = - (a^'^a^ + d^aa^a + a^a^a^ + dd^^d + dd^dd^ + d^d^^ 
W |d^d| = - (a^a + ac 



(63) 



6 
1 

2 

For the evolution operator we then flnd. 



(64) 
(65) 



U{t,to) =T7exp 



'to 



IK 
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2a){t')a{t') 



(66) 
(67) 



Relation (p2D also remains valid after replacing T+ — > and TCmtit) — > 7ij^(t). For any 
symmetrically ordered operator expression, W |-^}, we have 

(W {x}) = yrf2aiy(a)X|a(t)^„(t),at(t)^a*(t), 
where is the Wigner distribution characterising the initial state, 

Finally, we arrive at: 



(68) 



(69) 



(C-,C+,d,d) = Jd'aW{a)e^ce 



t I At, 



X e ^ T- T- 



I a_ (t)=a+ (t)=a(i),aL {t)=at {t)=a* (t) ' 



where 



(C-, C+, Cl Ci) = (t^ exp (C-at + CU) Tf exp (C+at + C^a 



(70) 



(71) 



Note that averages generated by functionals (0) and (^) may indeed differ. For instance, 
with h{t) = a){t)a{t), 



(r_n(t)T+n(t)) =([hit)] 



T^h{t) T^h{t) 



To obtain ([h{t)]^\ directly from ([ri|), we may use, e.g., the property 



[h{t)Y) = jim (T^a^(t)a(t + 6t) Tf a^(t + 6t)a{t) 
This limit should be preceded by the one associated with causal regularisation. 



(72) 
(73) 

(74) 



III. CAUSAL VARIABLES 



Our next goal is to investigate the causal structure of Eqs. ( |58D and ([70|). The concept of 
causality is introduced via the retarded Green's function, G{t). Following this, we are able 
to define the input and output of a quantum system. We then show that, physically, the 
input and output thus introduced correspond to a generalisation of Kubo's linear reaction 



approach [14] to a full nonlinear quantum-stochastic response problem. 



A. The method of normal ordering 



Consider in more detail the differential quadratic form in (pSD- Making use of relations 
(PH), and utilising the notation introduced by Eq. (|i3|), we find 

We now change the functional variables, a±(t),a±(t) — > a{t) , a'^ {t) , ^{t) , ^\t) , in order to 
obtain 



That is, 



' ^At> (78) 



6 _ 6 6 

6a{t) ~ 6a4t) ^ 6a.{t)' ^ ^ 

6 6 6 

+ ITTTTV (80) 



6a^it) 6al{t) 5aL(t)' 
These relations determine the new variables up to given functions which we chose to be zero: 

a+{t) = a{t), (81) 

a-{t) = a{t) -i^{t), (82) 

ai{t) = a^t)+ze{t), (83) 

al{t) = a\t). (84) 



Consider now the second exponent in (^8|): 

C-al + Cla_ + C+ai + Cla+ = (C- + C+) «^ + (Cl + Ct) a - tCU + (85) 
This clearly suggests another substitution, this time in the functional H itself: 

2(c-,c+,ci,cl) = <^(c,c^^,^^), (86) 
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where 



+zC+(t) = cr{t), 



(87) 
(88) 
(89) 
(90) 



and 



C+(t) = -tait), 

C-(t) =C(t)+^Or(t), 

Cl(t)=^at(t). 



(91) 
(92) 
(93) 
(94) 



In causal variables thus introduced, Eq. (|58| ) becomes: 

$ (C, C^ ^, = / d'aP{a) exp (^G^ + ^G*-^^ exp ((at + (^a + ^a^ + ^V) 

X exp S"!!!! a, a^) ,jt(t)=Q,*(t),^(t)=^t(t)=o , (95) 



where 



^int(e,e 



^, a, 



(96) 



It should be stressed that equation (|95D is universal, while all details of the problem enter 
through S'int- In general S'int is found as 



•^int a, a^^=i j dt \h (a\t), a{t) - i^{t)^ - h (a^(t) + i^\t), a{t) 



(97) 



where h is the normal representation of the interaction Hamiltonian, 



Hint = ■ h (a\ aj : . 



Generalisation of (P3| ) to multimode problems (see Sec. |V]) is also straightforward. 



(98) 
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B. The method of Wigner ordering 



Causal variables for the case of symmetric ordering are defined in a similar way. We 
require that 



6^ I^^I^ 5^' 



and then proceed as above to yield two sets of substitutions: 



and 



a+{t) = a{t) + -at), 
a_{t) = ait)-'-at), 

aUt) = a\t)-'-eit), 



Ciit) 



Ct(t)-^at(t), 



(99) 



(100) 
(101) 
(102) 
(103) 



(104) 
(105) 
(106) 
(107) 



(108) 



Further, on rewriting functional ( |7l| ) in the variables ( [104| )-( p!07D , 
we find that 

<l>^ (C, C^ cr, at) = J £aW{a) exp (^^^ + ^G*^^ exp ((at + + ia^ + ^a) 

X exp S*;^ (^^, ■, a, at) |„(j)^Q(^)^^t(t)=a*(t),5(t)=5t(t)=o- (109) 



Here, 



dt 



-h"^ (a\t) + -e{t),a{t) + -at) 



(110) 
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where is the symmetric representation of the interaction Hamiltonian, 



Hi„t = w{/i^ (a^a)}. 



;iii) 



For the Kerr nonhnearity, 



K 



t2^ 



(112) 



It should not be overlooked that the definition of causal variables is ordering-specific. In 
fact, the notation (, (\ a, a"^ is used for two different sets of causal variables: (|9lD -(p4D in 
the case of normal ordering, and (|104|) -( p!07|) in the case of symmetric ordering. Since it is 
always clear by the context which case is being discussed, no confusion should occur. 



C. Quantum nonlinear-reaction problem 

To gain more insight into the causal variables, we now introduce a variable pump term 
into the interaction Hamiltonian, which then reads: 

^int(t) = '^a^\t)a\t) + s(t)at(t) + s*{t)a{t). (113) 

The external source, s{t), is a given c-number function. We mark by tilde all quantitities 
defined in the presence of the source. With the source, 

^int = 5i„t + es* + e^s. (114) 



Moving ^s* + ^'s to the second exponent in (R^) results in an identity, 



$ (C, ^, cr^) = ^ (C, a + s,a^ + s*) , (115) 

which also holds for the <l>^/ <l>^ pair. 

This way, the physical information contained in the dependence of both $ (^(, (\ a, a^^ 
and $^ (^(, a, cr^^ on a and is the system's reaction to an external perturbations, while 
the variables a, o"^ define an input of the system. The fact that this input is also determined 
dynamically makes it independent of ordering. The variables (, define an output of the 
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system. Unlike the input, the output turns out to be ordering-specific. Consider firstly the 
case of normal ordering; assuming the source to be arbitrary, the full physical information 
about the system is obtainable from 

^{CXKo,o) = ^{CXKs,s*). (116) 

In turn, making use of (|5^) we get: 

I- (C, C^ 0, O) = exp (C^l) T+ exp ((a ) ) . (117) 

This is nothing but a characteristic functional of Glauber's renowned time-normal averages 
of the Heisenberg field operator in the presence of the source. The source terms in the 
Hamiltonian are also quite recognisable; they appear in Kubo's linear reaction theory | ]14|| . 
Introducing causal variables is thus equivalent to a nonlinear-reaction reformulation of a 
quantum system. 

Unlike the linear reaction theory, the nonlinear reaction theory depends explicitly on 
which quantities are to be "measured". Causal variables for the normal ordering, Eqs. 



(pT|)-(|9^), correspond to "measuring" time-normal averages of the field operators. Causal 
variables for the symmetric ordering, Eqs. ( |104|) - (|107| ), introduce another set of "measured" 
quantities. For example, relation ( |116| ) holds equally for $^ and $^ and in place of ( |117D , 



one finds 



(C, C^ 0, O) = (t^ exp (ic^a + l(i') Tf exp (^C^i + ^C^') ) • (118) 



Disregarding the time orderings turns ( p.l8| ) into a characteristic functional of symmetrically 
ordered operator averages. We shall therefore term the operator ordering introduced by 
( |118| ) a time-Wigner (Tw) ordering. That is, by definition, 

<l^(C,C^0,0) =^TH.exp(c^UCa)), (119) 

where 

T^exp (C^a + Ca^) = T_^exp gc^a+ ^C^^) Tf exp ^^(^a+^Ca^^ . (120) 
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We have dropped the tildes here to emphasise that, in itself, the definition of the time- 
Wigner ordering does not depend on the presence or absence of the sources (nor on any 
other details of the dynamics). 

Under time-normal ordering, the "most recent" creation (annihilation) operator becomes 
the leftmost (rightmost) in the product. The Tp^z-ordering acts in a similar manner, only 
symmetrised in respect of the creation and annihilation operators: 

Tw {x{t)m ■ ■ ■ m} = \ [Tw m') ■ ■ ■ Kn}m 

+ x{t)Tw{m---Kt")]]- (121) 
Here, x, y, ■ ■ ■ , z stand for field operators (i.e., a or a^), and t should exceed all other time 



arguments in the product, t > t\ - ■ ■ ,t". Equation (|121|) is a recurrence relation, defining 
Tw for the case of time arguments which are all different. For products of two operators, 
time-Wigner ordering is merely symmetric ordering, 

TwHmt') = I mm + mm = i m^mu ^ (122) 

where [■■■]+ denotes an anticommutator. For higher-order products, time ordering becomes 
essential, so that, for three operators, assuming that t > t',t", 

T^mmm) = \m. [mrm)uu- (123) 

It is easy to see that the following rule holds for an arbitrary Tvi/-ordered product: write 
all possible permutations of the operators in a product, then retain only those where time 
arguments are arranged in a C-contour sequence (i.e., times increase then decrease). 

Consider now the behaviour of time-Wigner averages for equal time arguments. Assum- 



ing t' > t, from ( |123| ) we find, 

T^mm^t' + 0) - r^x(t)j(t')z(t' - 0) = i[x(t), im, m]]+, im 

where [■ ■ ■] is a commutator, and Twx{t)y{t')z{t' ±0) = lim^t^o Twx(t)y(t')z{t' ±6t) (i.e. the 
limit applies to the Tiy-ordered product as a whole). This means that, unlike time-normal 
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products, time-Wigner products are not continuous at coinciding times, nor do the limits of 
Tiy-ordered products at equal times coincide with symmetric products. For instance, 

(125) 
(126) 

Symmetric averages may only be recovered as combinations of such limits from different 
directions, e.g.. 



1 






a 


~ 2 




1 






a 


~ 3 





1 2 

Wa\t)a^{t) = -Twa\t + {})8?{t) + -Twa\t)a{t)a{t + 0). 
3 3 



(127) 



IV. CALCULATING QUANTUM AVERAGES AS CLASSICAL STOCHASTIC 

AVERAGES 

In this section, we show that there exist classical stochastic problems, such that Eqs. 
(^) and ( |70D may be interpreted as stochastic averages. In particular, Eq. ( |58| ) yields a 
generalisation of the positive-P representation (well known in quantum optics) to a quantum 
nonlinear response problem. 



A. Classical stochastic response problem 

Relations (^) and ( |109| ) give formal solutions to quantum response problems (formu- 
lated, respectively, in terms of time- normal and time-Wigner averages). It is instructive to 
compare these relations to a solution to a classical stochastic response problem. To this end, 
consider a c-number stochastic field a(t), which obeys an integral equation, 

a(t) = a{t) + j dt'G{t - t')o^tot(t'), (128) 

where a(t) is the in-field and atot(t) is the field source. For the purposes of this paragraph, 
the kernel G(t) is assumed to be regular and retarded, G(t) = 0,t < 0, and otherwise 
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arbitrary. We assume that the in-field, a{t), is also arbitrary. (This allows one to regularise 
G{t) without changing a{t).) The full field source atot(^) consists of two parts, 

atot{t)=(T{t)+a'{t). (129) 

The external source, o"(t), is regarded as given, while the random source, cr'{t), depends 
on the field. That is, the random source describes the field's effective self-action (which 
physically originates, e.g., in interaction with a medium). As a random quantity, c'lt) is 
fully characterised by a probability distribution, conditional on the field at the same time t: 

U{a'{t)\a{t)). (130) 

Resolving formally the self-action problem results in a probability distribution over the 
random source, cr'{t), conditional on the in-field, a{t), and the external source, cr{t). This is 
found by substituting Eq. (|128|) for a(t). 



n{a'it)\[a + Gia + a')]{t)). (131) 

Importantly, this expression does not contain a vicious cycle because of the assumed regular- 
and-retarded nature of G(t): a'{t) depends on cr'(t') only for t' < t. 

Consider now statistical properties of the field. With the self-action resolved, these are 
also conditional on the in-field, a{t), and the external source, cr(t). For the characteristic 
functional of multi-time stochastic field averages we find, (with ([t) being an arbitrary 
function) 



S (CI a, a) = exp (^J dtC{t)a{t)^ = eC["+G(-+-')] (132) 
= J D°°(T'e^[°+^("+"')]n(aV + G'[a + (T']). (133) 

The upper bar here denotes an averaging over the statistics of a', which is afterwards ex- 
plicitly rewritten as a trajectorial (functional) integral. (Here and in what follows, we will 
make extensive use of a brief notation so as to prevent our relations from growing bulky.) 
Then, firstly, we pull G{<7 + a') out of 11 (o"'| a + G[o" + o"']) by applying a shift operator: 
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n (a'l « + G[(y + a']) = e^^^'^+'^'^H a) . (134) 

Secondly, we pull all factors except 11 out of the functional integral, resulting in: (with ^{t) 
being another arbitrary function) 

S(C|a,(T) =yD~(T'e^"+(^+^)^("+"')n((T'|a) (135) 
= e^"+('^+^)^K*) y'D-a'e«'^'n(a'|a)|^=o (136) 
= e^"+(«+^)^K*)e^(^l")|^=o. (137) 

We have introduced a characteristic functional of cumulants of the random source conditional 
on the full field, 

S{C\a) = lnjD°°a'e^'''U{a'\a). (138) 
After further algebra, which is much assisted by Eq. (^51), we arrive at 

E(C|«,cr) =e^^*e^'^+«'^e^(^l'^)|5=o,a=a. (139) 

Disregarding the averaging over the pseudo-distributions, Eqs. ( P^ ) and ( |1U9D look very 
much like generalisations of Eq. (|139|) to a pair of random fields. It is not however immedi- 
ately obvious that S'int and may be interpreted as characteristic functional of cumulants 
of some classical noises. In fact, this turns out to be unconditionally the case for S'int, whereas 
for this interpretation holds only with discretisation of the time axis. 

B. Hubbard-Stratonovich transformation: introducing noise sources constructively 

1. Generalising the positive-P representation to a quantum problem of nonlinear reaction 

Consider a standardised real delta-correlated Gaussian noise, x(t)y 

X{t)xit')=5it-t'). (140) 

The following relation (Hubbard-Stratonovich Transformation — HST) holds for an arbitrary 
function, x(t), 
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so we find, 



exp S'int C\ CL, a^) = exp 



(142) 



where xK't) is another standardised real delta-correlated Gaussian noise, uncorrelated with 
x(t), and the upper bar denotes averaging over the statistics of the noises. The exponent on 
the RHS of ( |142| ) is linear in ^(t) and ^^(t). This means that for given x{t) and x''{t)y the 
equations for a{t) and a''"(t) become regular, with the random sources being, respectively, 



In other words, a{t) and a'^{t) obey a pair of coupled stochastic integral equations. 



ait) + '^a\t')a\t') + f-^xit')a{t') 



(143) 
(144) 



a{t) = a{t) + J dt'G{t - t') 
3t(t) =a*{t) + j dt'G*{t-t') 



a\t) + ^a^\t')a{t') + J^x\t')a\t') 



(145) 
(146) 



On removing the causal regularisation of G{t), one recovers a pair of Ito stochastic differential 
equations, 

. da{t) 



<y{t) + \a\t')a\t') + J'-^x{t')a{t'\ 



-IK 



(147) 
(148) 



dt ' '2 ' ''' ' ■ V 2 

The fact that Ito calculus should be chosen is due to the causal regularisation of G'(t), which 
makes sources at time t independent of fields at the same time, which is the characteristic 
property of Ito calculus. 

Without the external sources (i.e., with a = = 0), equations ( |147| ), ( |148D are well 
known in quantum optics under the name of the positive-P representation. In that form, 
they allow one to calculate time-normal averages of the field operators, describing radiation 
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properties of the system. With sources included, they become apphcable to a much wider 
assemblage of operator averages, covering the full quantum-stochastic nonlinear reaction 
problem. 



2. Positive- W representation 

At a first glance, factorising noise terms in S",^ does not pose a problem either. Intro- 
ducing a standardised complex delta-correlated Gaussian noise, ri{t), such that 



7]{t)7^{t') = 0, 7^{t)7^*{t') = 6{t - t'), (149) 

allows one to factorise arbitrary products: 



Qxy^^^n+yr_ (150) 
We shall write the real and complex HST's, Eqs. (|141|) and ( p.5CI|) , as: 



2 

XX, xy =^ xrj + yrj . (l^lj 



2 

These relations imply the definition of the corresponding noises as standardised Gaussian 
noises (respectively, real and complex ones). Then, e.g., (cf. Eq. ( |112| )) 

- ^e^a ^V^-^- V'^^a ^ xVv^' " V*1^a. (152) 



Consider now the random quantity, z/(t) = xi^)yvi^) ■ x(^) and ri{t) are independent, and 
we have. 



iy{t)u*{t') = x{t)x{t') ^/v{t)v*{t') = S{t - t')^7r6{t-t') = 6{t - (153) 

That is, we have encounted an infinity (divergence). The underlying reason is certainly 
that both x{t) and ri{t) are highly singular functions and the "quantity" z/(t) is simply not 
defined. 

For all practical purposes (such as computer simulations) stochastic differential equa- 
tions will be replaced by finite-difference equations. We therefore assume an equidistant 
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discretisation of the time axis, with At being the step of the time grid. The integration and 
the (5-function are understood as usual, 

fdt = AtJ2, 6{t-t')=At-'6u', (154) 

where Su' is the Kronecker symbol on the time grid (i.e. 6u' = 1 if t = t' and otherwise). 
We renormalise the standardised Gaussian sources, 

V{t) = At-'/'vit), x{t) = Ar'/'m, (155) 
making them 6tt' (Kronecker symbol) correlated, not 6{t — t') (5-function) correlated: 



rj{t)r]*{t')=6u', xmit')=6u'. (156) 

The short upper bar always marks the Kronecker-correlated noises (and should not be con- 
fused with the long one denoting the averaging; the latter always extends over the whole 
expression). The HSTs then become, 

x'^ X XX _^xf] + yf]* 

xy =^ -1= — . (157) 



'At VAt 

There are at least two ways in which the third order terms may be factorised. One of 
them, generalising ( |152|) , employs two complex and two real HST's: 

--^t2^a ^ pAtf]^ + qf]X xVP^i'^ + (158) 



2 

2 



where the functions p(t), q{t),p^{t), q^{t) obey the relations 

and are otherwise arbitrary (they can even be random, given they are not correlated with the 
Gaussian noises, fjit), x(t), f7"l'(t), xK't))- This results in the following coupled finite-difference 
equations for a{t) and a^{t): 
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—I 



i [a{t + At) - a{t)] 



a\t)a{t) 
a\t)a{t) 



a(t) At, 
al'(t) At, 



(161) 
(162) 



where 



/it(t) = xt(t)v5Ww+?wrw- 



(163) 
(164) 



If we require the average random excursion squared at each time step, 



At^ 



At' 



\p{t)\+ pHt)]+mf+ qHt) 



(165) 



to be minimal, we find. 



71 K \a{t)\ 
8At2 



8At2 



(166) 



Apart from Eq. (|16(]|)), phases of the p's and g's are of no physical consequence. We choose: 



q\t) = 
For ( p.65| ), we then find. 



TT K,a{t) 
8At2 

8At2 



1/3 



1/3 



Pit) 



pt(t) 



71 na(t) 



8At2 



2/3 



8At2 



2/3 



(167) 
(168) 



At^ 



IMt)l' + |/it(t)|^ 



oc 



At2/3. 



(169) 



This should be compared with the scaling characteristic of a Wiener process, oc At. The 
"shortage" of At^^^ results in an overall random excursion squared over time T scaling as, 
oc T/At^/'^. This is indeed better that the scaling, oc T/At^/^, characteristic of the "naive" 
factorisation, Eq. ( |152|) , yet the improvement is only quantitative, the process still diverging 
in the continuous time limit. From a practical perspective, this means that sampling noise 
increases with decreasing time step. 

Another way of factorising the third order terms employs three complex HSTs: 
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(170) 
(171) 
(172) 



As above, the p's and g's here are arbitrary functions obeying Eq. (|160D , and r(t)r^(t) = 1/2. 
That is, we recover the generic equations (|158|) and ( |159| ), this time with 



/i(t) = rrj'yjpr] + p^t]'^ + qrj*, 



(173) 

(174) 



(Note that the terms proportional to the gs have exchanged places.) Minimising the random 
excursion squared, Eq. ( |165| ), we find the optimised values of free parameters to be 



r(t) = r\t) = 



V2' 



64At4(|a(t)| + |at(t)|) 
7r/t2at3(t) 



1/6 



nl/6 



(175) 
(176) 
(177) 



_64At4(|a(t)| + |at(t)|)_ 

with the ps being recovered from Eq. ( |16CI| ). Evidently both ways of factorising the third 
order terms result in the same scaling of the noises. 



C. Numerical experiment 



We have calculated the time dependence of the quadrature amplitude. 



X{t) = (hit) + a) it)) = a{t) + at(t). 



(178) 



assuming that the initial state of the oscillator is a coherent state \a). An analytical expres- 
sion for this quantity may be found in [Q. In Fig. ^, we plot the relative error in X{t), 

^simulated ('^) ^exact (^) 



sx 



^cxact (^) 



(179) 
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where Xsimuiatcd(^) is found via stochastic simulations of the "+W" equations, ( |161|) - (|162|) , 



for K = 4 and a = 1. To have an idea how important the third order noises are, we also 
calculated X(t) using the truncated Wigner representation [0. The latter is obtained by 
dropping the noise sources in Eqs. ( |161| )-( p^ ). The resulting equations are readily solved 
analytically; the averaging over the initial W-distribution does not pose a problem. The 
relative error in X{t) found via the truncated Wigner ("— W") representation is also shown 
in Fig. ^. As is evident from inspection of Fig. |^, the "+W" result is a manifest improvement 
over the "— W" one. The bad news is, however, that the sampling error remains rather large, 
even after averaging over more than half a billion trajectories. A further piece of bad news is 
that simulations fell victim to numerical instability shortly after the maximal time shown in 
Fig. 1^. These problems notwithstanding, the very fact that third order noises my be subject 
to stochastic simulations has been clearly established. 

V. THE COOKBOOK 

A. The general recipe 

We now summarise our results, recipe-style, the way they should be applied to practical 
calculations. For an n-mode system, Schrodinger-picture annihilation and creation operators 
are vectors in respect of the mode index, 

as = {dsk}, a^ = {4fc}' k = l,---,n, (180) 

the interaction-picture field operators, 

a(t) = {afc(t)}, a\t) = {dlit)} , k = l,---,n. (181) 

The system Hamiltonian consists, as usual, of the free and interaction Hamiltonians, cf. Eq. 
(HI). The free Hamiltonian is an arbitrary matrix, H = {Hkk'}, k, k' = 1, ■ ■ ■ ,n, in the mode 
indices. The free Schrodinger equation thus reads 

.da{t) 



dt 
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Ha{t). (182) 



Quantum-classical mappings yield the generic system of 2 x n equations, 

= Ha{t) + (T{t) + (T'{t), (183) 



184) 



for 2 X n c-number random fields, 

Kt) = {afc(t)}, at(t) = {4(t)}, A; = l,---,n. (185) 



The vectors of random sources, cr'(t), cr'"^(t), depend on the interaction and on the actual 
type of operator ordering underlying the mapping. 

Consider firstly the case of time- normal ordering. Without the given sources (i.e., with 
cr{t) = cr'i'(t) = 0), stochastic averages of the random fields coincide with the quantum 
averages of the time-normally-ordered products of the Heisenberg field operators, a.k(t), a^(t), 
k = 1, ■ ■ ■ ,n. With the given sources, Eqs. ( |183| )-( [TS^ ) constitute a constructive nonlinear- 



reaction formulation of the quantum system. The random sources are deduced from the 
normal form of the interaction Hamiltonian, 

TYint = ■.h[as,al) : . (186) 

Namely, one should calculate the functional, (cf. Eq. (|96|)) 

5int |^ a, a^) = ij dt [h (at(t), a{t) - i^{t)) - h [a\t) + t^\t), a{t))] . (187) 

Since Tiint is Hermitian, this can be further simplified resulting in 

Sint ^^ a, a^) = ijdth (a^(t), a{t) - z^(t)) + conj., 

where conjugation acts as a formal Hermitian transformation (i.e., it interchanges quantities 
with and without dagger, a{t) ^ a^it), ^{t) ^ ^^(^), and complex-conjugates other c- 
numbers). One should then factorise powers and products of ^s, employing a suitable set of 
Hubbard-Stratonovich transformations (see Eqs. ( |151D ), until an expression emerges which 
is already linear in respect of the ^s: 
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n 



k=l 



(188) 



(Note that terms without in ^int always cancel) The as thus obtained are exactly those 
which appear in Eqs. (|183D -( p!8^) . If S'int is at most quadratic in ^'s, Eqs. ( [183|) - (|184|) 
are a system of genuine Ito stochastic differential equations. Otherwise they can only be 
interpreted as difference equations over a finite time step At. 

In the case of the time-Wigner ordering, this recipe applies with the replacement S'int — ^ 
S^, where 



SZ a, a'^)=ij dt [h^^ (at(t) - '-^\t), a(t) - ^^(t) 



h^{a^it)+'-eit),ait)+'-m 



\jdth^ (a^t) - '-i\t), a{t) - '-i{t)^ + conj., 



(189) 
(190) 



and h}^ is the symmetric form of the interaction Hamiltonian, 



?T:int = w{/i^ (as,4)}- 



;i9i) 



We are however not aware of any nonlinear system where SZ would turn out to be quadratic 



in the ^'s. 



B. Degenerate OPO 

1. Generalised positive-P representation 

We shall now illustrate the general recipe by deriving the positive-P and positive-W 
representations for a degenerate optical parametric oscillator (OPO). The degenerate OPO 
consists of two coupled oscillators described by the Hamiltonian, 

n = u;asi4i + 2u;as2a^2 + y «si«S2 - ali«s2 • (192) 
Then, in brief notation. 
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af (aa - ^6) - (^i - ^^lY 4 



+ conj. 



-m^iai4 + ^^64^ - -^^14 + conj. 



— i^i |^/taia2 + X V 'f^" 
+<2-ai , 



(193) 



where X^(^) ^ire a pair of independent real 5-correlated Gaussian noises. The positive-P 
representation for the OPO, generahsed to the response problem, then reads, 



da\{t) 

dt 

da2{t) 



dt 
dalit) 
dt 



-ia2{t) - -4(t), 



(194) 

(195) 
(196) 
(197) 



This derivation is strikingly simple and straightforward, and compares very favorably to the 
common derivation based on phase-space techniques (not to mention that Eqs. ( |194|) - 



|197|) allow for a much deeper insight). 



2. Positive- W representation 



Similarly, 



SZ ^^ a, at) = Ua\ - ^^1/2)' (as - ^6/2) - (ai - ^6/2)' (4 - ^^1/2 



.IK 

-z<iaia^ + ^-4201 + -^?i^2 + conj. 

! ■ t t I ■'^f t2 *f2 , ■ ft , 

-z<iaia^ + z-42ai - + ^^'^G + conj. 



it* 



-iix ( fi:ai4 + X^\/pV 



+ conj. 
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K 



K 



resulting in the positive-W representation for the OPO: 

= -^ai(t) + Kal{t)a2{t) + x{t) p{tW* (t) , 



dt 
da2{t) 

dt 

dalit) 

dt 



K 



-icy2{t)-'-al{t) + q{t)r^{t), 
Mal{t)-^a\\t) + q\tW{t). 



(198) 

(199) 
(200) 
(201) 
(202) 



Here, rj{t),rii\t) and xif)iX^if) are pairs of standardised Gaussian noises (respectively, com- 
plex and real ones), and the functions p(t),p^(t), g(t), q\t) obey the relations. 



p{t)q\t) = p\t)q(t) 



(203) 



Equations ( |199| )- (|20^) are difference equations over a finite time step dt = At, with the 
noises normalised such that 

5tt' 



X{t)xit') = xKt)xKt') = vitWit') = r]^{t)r]^*{f) = Sit - 1') 



(204) 



The p's and g's will be chosen so as to minimise the weighted sum of random increments 
squared. 



de 



Xit)^p{t)r,^*{t) + xKt)^PKtWit) 



+ w 



p{t)V7rdi + 



WK^dt 



\q{t)m" + \qKt)vKt)? 



conj. 



(205) 



16p2(t) 

where we have chosen the p's to be real and positive. The weighting factor u; > allows 
one to redistribute the noise between the fundamental and the subharmonic, and may be 
fine-tuned by trial and error when running stochastic simulations. Minimising (|205|) yields: 



p{t)=p\t) = ]-^WK''^dt/^ 



(206) 



The g's are then recovered from (|203|) . Accounting for (|206|) , all noise increments in the 
positive-W equations scale as, cx dt^^'^. 
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FIGURES 
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FIG. 1. The Schwinger-Perel-Keldysh C-contour (thin hnes) and the three contractions (thick 
Unes) contributing to Ac, cf. Eq. (^). The arrows on contractions are from creation to annihilation, 
cf. Eq. (p7[), and the time order of ends corresponds to contractions being non-zero. 
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FIG. 2. Relative error in the coherent field amplitude X{t), calculated using the positive-W 
(+W, solid line) and truncated Wigner (— W, dash-dotted line) representations, for nonlinearity 
K = 4, and the inital coherent state of the oscillator with amplitude a = 1. 
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